#include "wrapper_hypre.h"
#include <vector>

typedef HYPRE_Int idx_t;
typedef HYPRE_Real data_t;

const idx_t glb_nparts = 4;
const idx_t num_diag = 19;

int main(int argc, char* argv[])
{
    int my_pid, num_procs;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

    idx_t arg_cnt = 1;
    const idx_t glb_dims [3] = { atoi(argv[arg_cnt++]), atoi(argv[arg_cnt++]), atoi(argv[arg_cnt++]) };
    const idx_t prc_dims [3] = { atoi(argv[arg_cnt++]), atoi(argv[arg_cnt++]), atoi(argv[arg_cnt++]) };// 注意这是块内的进程数
    assert(glb_dims[1] % glb_nparts == 0);
    const idx_t blk_dims[3] = { glb_dims[0], glb_dims[1] / glb_nparts, glb_dims[2] };

    const std::string pathname(argv[arg_cnt++]);
    const std::string solver_name = std::string(argv[arg_cnt++]),
                        prec_name = std::string(argv[arg_cnt++]);

    const int nprocs_per_part = prc_dims[0] * prc_dims[1] * prc_dims[2];
    const idx_t my_part = my_pid / nprocs_per_part,
                pid_in_part = my_pid - my_part * nprocs_per_part;
    idx_t cart_ids[3];
    {// 计算本进程在part内处于什么位置
        cart_ids[0] =  pid_in_part / (prc_dims[1] * prc_dims[2]);
        cart_ids[1] = (pid_in_part - cart_ids[0] * prc_dims[1] * prc_dims[2]) / prc_dims[2];
        cart_ids[2] = (pid_in_part - cart_ids[0] * prc_dims[1] * prc_dims[2] - cart_ids[1] * prc_dims[2]);
    }

    idx_t my_ilower[3], my_iupper[3], box_ends[3], my_dims[3];
    idx_t my_nelems = 1;
    for (idx_t d = 0; d < 3; d++) {
        my_ilower[d] =  cart_ids[d]     * (blk_dims[d] / prc_dims[d]);
        my_iupper[d] = (cart_ids[d] + 1)* (blk_dims[d] / prc_dims[d]) - 1;
        box_ends[d] = my_iupper[d] + 1;
        my_dims[d] = my_iupper[d] - my_ilower[d] + 1;
        my_nelems *= my_dims[d];
    }
#ifdef _zyDEBUG
    // for (int p = 0; p < num_procs; p++) {
    //     if (my_pid == p) {
    //         printf("Proc %d part %d: (%d,%d,%d)~(%d,%d,%d)\n", my_pid, my_part,
    //             my_ilower[0], my_ilower[1], my_ilower[2],
    //             my_iupper[0], my_iupper[1], my_iupper[2] );
    //         fflush(stdout);
    //     }
    //     MPI_Barrier(MPI_COMM_WORLD);
    // }
    // MPI_Barrier(MPI_COMM_WORLD);
#endif

    // 先读入本进程所负责的区域的数据
    data_t * A_buf = new data_t [my_nelems * num_diag];
    data_t * rhs_buf = new data_t [my_nelems];// 读入数据时从外到内为 [3360][3072]
    data_t * x0_buf  = new data_t [my_nelems];// 读入数据时从外到内为 [3360][3072]

    {// 读入
        float * vec_buf = new float [my_nelems];
        data_t my_x_dot = 0.0, my_b_dot = 0.0;

        MPI_File fh = MPI_FILE_NULL;// 文件句柄
        const MPI_Datatype etype = MPI_FLOAT;
        MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
        idx_t size[3], subsize[3], start[3];
        size   [0] = glb_dims [0]; size   [1] = glb_dims [1]; size   [2] = glb_dims [2];
        subsize[0] = my_dims  [0]; subsize[1] =  my_dims [1]; subsize[2] =  my_dims [2];
        start  [0] = my_ilower[0]; start  [1] = my_ilower[1] + my_part * glb_dims[1] / glb_nparts;// 注意这里
        start  [2] = my_ilower[2];

        MPI_Type_create_subarray(3, size, subsize, start, MPI_ORDER_C, etype, &read_type);
        MPI_Type_commit(&read_type);
        MPI_Offset displacement = 0;
        displacement *= sizeof(*vec_buf);// 位移要以字节为单位！
        MPI_Status status;

        int ret;
        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_b").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open b file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, vec_buf, my_nelems, etype, &status);
        MPI_File_close(&fh);
        #pragma omp parallel for collapse(2) schedule(static) reduction(+:my_b_dot)
        for (idx_t i0 = 0; i0 < my_dims[0]; i0++)
        for (idx_t i1 = 0; i1 < my_dims[1]; i1++)
        for (idx_t i2 = 0; i2 < my_dims[2]; i2++) {
            data_t tmp = vec_buf[((i0 * my_dims[1]) + i1) * my_dims[2] + i2];
            rhs_buf[((i2 * my_dims[1]) + i1) * my_dims[0] + i0] = tmp;
            my_b_dot += tmp * tmp;
        }
        // #pragma omp parallel for schedule(static)
        // for (idx_t e = 0; e < my_nelems; e++)
        //     rhs_buf[e] = vec_buf[e];

        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_x").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open x file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, vec_buf, my_nelems, etype, &status);
        MPI_File_close(&fh);
        #pragma omp parallel for collapse(2) schedule(static) reduction(+:my_x_dot)
        for (idx_t i0 = 0; i0 < my_dims[0]; i0++)
        for (idx_t i1 = 0; i1 < my_dims[1]; i1++)
        for (idx_t i2 = 0; i2 < my_dims[2]; i2++) {
            data_t tmp = vec_buf[((i0 * my_dims[1]) + i1) * my_dims[2] + i2];
            x0_buf[((i2 * my_dims[1]) + i1) * my_dims[0] + i0] = tmp;
            my_x_dot += tmp * tmp;
        }
        // #pragma omp parallel for schedule(static)
        // for (idx_t e = 0; e < my_nelems; e++)
        //     x0_buf[e] = vec_buf[e];
        
        // 读入SOA存储的矩阵
        for (idx_t d = 0; d < num_diag; d++) {// 逐条对角线
            ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_a." + std::to_string(d)).c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
            if (ret != MPI_SUCCESS) {
                printf("Could not open x file: ret %d\n", ret);
            }
            MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
            MPI_File_read_all(fh, vec_buf, my_nelems, etype, &status);
            MPI_File_close(&fh);
            #pragma omp parallel for collapse(2) schedule(static)
            for (idx_t i0 = 0; i0 < my_dims[0]; i0++)
            for (idx_t i1 = 0; i1 < my_dims[1]; i1++)
            for (idx_t i2 = 0; i2 < my_dims[2]; i2++) {
                idx_t src_e = ((i0 * my_dims[1]) + i1) * my_dims[2] + i2;
                idx_t dst_e = ((i2 * my_dims[1]) + i1) * my_dims[0] + i0;
                A_buf[dst_e * num_diag + d] = vec_buf[src_e];
            }
            // #pragma omp parallel for schedule(static)
            // for (idx_t e = 0; e < my_nelems; e++)
            //     A_buf[e * num_diag + d] = vec_buf[e];
        }
        delete [] vec_buf;

        MPI_Type_free(&read_type);
        data_t glb_b_dot = 0.0, glb_x_dot = 0.0;
        MPI_Allreduce(& my_b_dot, & glb_b_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(& my_x_dot, & glb_x_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        if (my_pid == 0) {
            printf("(myb, myb) = %.20e\n",  (double)glb_b_dot);
            printf("(myx, myx) = %.20e\n",  (double)glb_x_dot);
        }
    }

    idx_t obj_type = HYPRE_SSTRUCT;
    if (prec_name == "AMG") obj_type = HYPRE_PARCSR;
    if (my_pid == 0) {
        printf("object type : %d\n", obj_type);
    }

    HYPRE_Init();

    HYPRE_SStructGrid ssgrid;
    HYPRE_SStructStencil stencil;
    HYPRE_SStructGraph ssgraph;
    HYPRE_SStructMatrix   A;
    HYPRE_SStructVector   b, x, y, x0;
    HYPRE_ParCSRMatrix par_A;
    HYPRE_ParVector par_b, par_x, par_y;
    TEST_RECORD records[TEST_CNT];
    TEST_CONFIG config;
    config.case_name = "gfs";
    config.atol = 1.0e-6;
    
    {
        idx_t offsets[num_diag][3] = {
            // { 0, -1, -1},
            // {-1,  0, -1}, { 0, 0, -1}, { 1, 0, -1},
            // { 0,  1, -1},

            // {-1, -1, 0}, { 0,-1, 0}, { 1,-1, 0},
            // {-1,  0, 0}, { 0, 0, 0}, { 1, 0, 0},
            // {-1,  1, 0}, { 0, 1, 0}, { 1, 1, 0},

            // {  0,-1, 1},
            // { -1, 0, 1}, { 0, 0, 1}, { 1, 0, 1},
            // {  0, 1, 1},

            {-1,-1, 0},
            {-1, 0,-1}, {-1, 0, 0}, {-1, 0, 1},
            {-1, 1, 0},

            { 0,-1,-1}, { 0,-1, 0}, { 0,-1, 1},
            { 0, 0,-1}, { 0, 0, 0}, { 0, 0, 1},
            { 0, 1,-1}, { 0, 1, 0}, { 0, 1, 1},

            { 1,-1, 0},
            { 1, 0,-1}, { 1, 0, 0}, { 1, 0, 1},
            { 1, 1, 0},
        };
        struct COO {
            idx_t my_index[3];
            idx_t ngb_bid = -1;
            idx_t ngb_index[3];
            data_t val = 0.0;
        };
        std::vector<COO> offd_nnzs;
        for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
        for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
        for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
            data_t* dst_ptr = A_buf + (((i2 - my_ilower[2] ) * (my_iupper[1] - my_ilower[1] + 1)
                                    +    i1 - my_ilower[1] ) * (my_iupper[0] - my_ilower[0] + 1)
                                    +    i0 - my_ilower[0] ) * num_diag;

            COO possi;
            possi.my_index[0] = i0;
            possi.my_index[1] = i1;
            possi.my_index[2] = i2;
            
            for (idx_t d = 0; d < num_diag; d++) {
                const idx_t * ngb_offsets = offsets[d];
                idx_t ngb_i[3] = {  i0 + ngb_offsets[0],
                                    i1 + ngb_offsets[1],
                                    i2 + ngb_offsets[2] };
                idx_t ngb_part = -1;
                idx_t to_id[3];

                if (0 <= ngb_i[0] && ngb_i[0] < blk_dims[0] &&
                    0 <= ngb_i[1] && ngb_i[1] < blk_dims[1] ) {
                    if (ngb_i[2] < 0 || ngb_i[2] >= blk_dims[2]) {
                        assert(dst_ptr[d] == 0.0);
                    }
                    continue;
                } else {
                    if (ngb_i[2] < 0 || ngb_i[2] >= blk_dims[2]) {
                        assert(dst_ptr[d] == 0.0);
                        continue;
                    }

                    if (0 <= ngb_i[0] && ngb_i[0] < blk_dims[0]) {// 无需跨极点，普通的周期边界
                        // 则必然是ngb_i[1]超出范围了
                        if (ngb_i[1] < 0) { assert(ngb_i[1] == -1);
                            ngb_part = (my_part - 1 >= 0) ? (my_part - 1) : glb_nparts - 1;
                            ngb_i[1] = blk_dims[1] - 1;
                        } else { assert(ngb_i[1] == blk_dims[1]);
                            ngb_part = (my_part + 1 < glb_nparts) ? (my_part + 1) : 0;
                            ngb_i[1] = 0;
                        }
                    } else {// 跨极点
                        idx_t oppo_part = (my_part + glb_nparts / 2) % glb_nparts;// 正对面的块
                        if (ngb_i[0] < 0) { assert(ngb_i[0] == -1);
                            ngb_i[0] = 0;
                        } else { assert(ngb_i[0] == blk_dims[0]);
                            ngb_i[0] = blk_dims[0] - 1;
                        }
                        if (ngb_i[1] < 0) { assert(ngb_i[1] == -1);
                            ngb_part = (oppo_part - 1 >= 0) ? (oppo_part - 1) : glb_nparts - 1;
                            ngb_i[1] = blk_dims[1] - 1;
                        } else if (ngb_i[1] >= blk_dims[1]) { assert(ngb_i[1] == blk_dims[1]);
                            ngb_part = (oppo_part + 1 < glb_nparts) ? (oppo_part + 1) : 0;
                            ngb_i[1] = 0;
                        } else {
                            ngb_part = oppo_part;
                        }
                    }

                    possi.ngb_bid = ngb_part;
                    possi.ngb_index[0] = ngb_i[0]; assert(ngb_i[0] >= 0 && ngb_i[0] < blk_dims[0]);
                    possi.ngb_index[1] = ngb_i[1]; assert(ngb_i[1] >= 0 && ngb_i[1] < blk_dims[1]);
                    possi.ngb_index[2] = ngb_i[2]; assert(ngb_i[2] >= 0 && ngb_i[2] < blk_dims[2]);
                    possi.val = dst_ptr[d];
                    dst_ptr[d] = 0.0;// 清零

                    offd_nnzs.emplace_back(possi);
                }
            }// d loop
        }// elem loop
        const idx_t loc_offd_nnz = offd_nnzs.size();
#ifdef DEBUG_PRINT
        for (int p = 0; p < num_procs; p++) {
            if (my_pid == p) {
                printf("Proc %d offd nnz %d\n", my_pid, loc_offd_nnz);
                fflush(stdout);
            }
            MPI_Barrier(MPI_COMM_WORLD);
        }
        MPI_Barrier(MPI_COMM_WORLD);
#endif
        idx_t glb_offd_nnz = 0;
        MPI_Reduce(& loc_offd_nnz, & glb_offd_nnz, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
        if (my_pid == 0) {
            printf("Global Offd nnz %d\n", glb_offd_nnz);
        }

        HYPRE_SStructGridCreate(MPI_COMM_WORLD, 3, glb_nparts, &ssgrid);
        HYPRE_SStructGridSetExtents(ssgrid, my_part, my_ilower, my_iupper);

        HYPRE_SStructVariable vtypes[1] = { HYPRE_SSTRUCT_VARIABLE_CELL };
        for (idx_t part = 0; part < glb_nparts; part++) {// 每个part都要执行定义变量的操作
            HYPRE_SStructGridSetVariables(ssgrid, part, sizeof(vtypes)/sizeof(HYPRE_SStructVariable), vtypes);
        }
        HYPRE_SStructGridAssemble(ssgrid);// a collective call finalizing the grid assembly.
        
        HYPRE_SStructStencilCreate(3, num_diag, &stencil);
        for (idx_t e = 0; e < num_diag; e++) {
            HYPRE_SStructStencilSetEntry(stencil, e, offsets[e], 0);// 只有0号变量一种
        }

        HYPRE_SStructGraphCreate(MPI_COMM_WORLD, ssgrid, &ssgraph);
        HYPRE_SStructGraphSetObjectType(ssgraph, obj_type);
        for (idx_t part = 0; part < glb_nparts; part++) {// 每个part都要执行定义模板的操作
            HYPRE_SStructGraphSetStencil(ssgraph, part, 0, stencil);
        }

        // 添加上非stencil的连接关系（跨极点连接关系）
        for (size_t k = 0; k < offd_nnzs.size(); k++) {
            COO & possi = offd_nnzs[k];
            HYPRE_SStructGraphAddEntries(ssgraph, my_part, possi.my_index, 0, possi.ngb_bid, possi.ngb_index, 0);
        }
        HYPRE_SStructGraphAssemble(ssgraph);

        // 建立向量
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &b);// Create an empty vector object
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &y);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x0);
        HYPRE_SStructVectorSetObjectType(b, obj_type);// Set the object type for the vectors to be the same as was already set for the matrix
        HYPRE_SStructVectorSetObjectType(x, obj_type);
        HYPRE_SStructVectorSetObjectType(y, obj_type);
        HYPRE_SStructVectorSetObjectType(x0, obj_type);
        HYPRE_SStructVectorInitialize(b);// Indicate that the vector coefficients are ready to be set
        HYPRE_SStructVectorInitialize(x);
        HYPRE_SStructVectorInitialize(y);
        HYPRE_SStructVectorInitialize(x0);
        {// 填充向量数据
            HYPRE_SStructVectorSetBoxValues(b, my_part, my_ilower, my_iupper, 0, rhs_buf);
            HYPRE_SStructVectorSetBoxValues(x, my_part, my_ilower, my_iupper, 0, x0_buf);
            HYPRE_SStructVectorSetConstantValues(y, 0.0);
        }
        HYPRE_SStructVectorAssemble(b);
        HYPRE_SStructVectorAssemble(x);
        HYPRE_SStructVectorAssemble(y);// 对y也要创建
        HYPRE_SStructVectorAssemble(x0);
        HYPRE_SStructVectorCopy(x, x0);

        HYPRE_SStructMatrixCreate(MPI_COMM_WORLD, ssgraph, &A);// Create an empty matrix object
        HYPRE_SStructMatrixSetObjectType(A, obj_type);
        HYPRE_SStructMatrixInitialize(A);// Indicate that the matrix coefficients are ready to be set 
        idx_t stencil_indices[num_diag];
        for (idx_t j = 0; j < num_diag; j++)
            stencil_indices[j] = j;
        {// 填充矩阵数据
            // 结构部分
            HYPRE_SStructMatrixSetBoxValues(A, my_part, my_ilower, my_iupper, 0, num_diag, stencil_indices, A_buf);
            idx_t * nentries = new idx_t [my_nelems];// 记录已经为该位置设置了多少个非零元 
            #pragma omp parallel for schedule(static)
            for (idx_t i = 0; i < my_nelems; i++) nentries[i] = num_diag;
            // 非结构部分
            for (size_t k = 0; k < offd_nnzs.size(); k++) {
                idx_t * my_index = offd_nnzs[k].my_index;
                const idx_t loc_1D = ((my_index[2] - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                    +  my_index[1] - my_ilower[1]) * (my_iupper[0] - my_ilower[0] + 1)
                                    +  my_index[0] - my_ilower[0];
                HYPRE_SStructMatrixSetValues(A, my_part, my_index, 0, 1, & nentries[loc_1D], & offd_nnzs[k].val);
                nentries[loc_1D] ++;// 计数增加
            }
        }// fill in matrix data
        HYPRE_SStructMatrixAssemble(A);// a collective call finalizing the matrix assembly.


        for (idx_t test = 0; test < TEST_CNT; test++) {
            HYPRE_SStructVectorCopy(x0, x);
            
            data_t b_dot, x_dot, Ab_dot, Ax_dot;
            idx_t ret;        
            // 做spmv检验一下数据传对没有
            if (obj_type == HYPRE_PARCSR) {
                HYPRE_SStructMatrixGetObject(A, (void **) &par_A);
                HYPRE_SStructVectorGetObject(b, (void **) &par_b);
                HYPRE_SStructVectorGetObject(x, (void **) &par_x);
                HYPRE_SStructVectorGetObject(y, (void **) &par_y);

                ret = HYPRE_ParVectorInnerProd(par_b, par_b, &b_dot);
                ret = HYPRE_ParVectorInnerProd(par_x, par_x, &x_dot);
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_b, 0.0, par_y);
                ret = HYPRE_ParVectorInnerProd(par_y, par_y, &Ab_dot);
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_x, 0.0, par_y);
                ret = HYPRE_ParVectorInnerProd(par_y, par_y, &Ax_dot);
            } else { assert(obj_type == HYPRE_SSTRUCT);
                ret = HYPRE_SStructInnerProd(b, b, &b_dot);
                ret = HYPRE_SStructInnerProd(x, x, &x_dot);
                HYPRE_SStructMatrixMatvec(1.0, A, b, 0.0, y);
                ret = HYPRE_SStructInnerProd(y, y, &Ab_dot);
                HYPRE_SStructMatrixMatvec(1.0, A, x, 0.0, y);
                ret = HYPRE_SStructInnerProd(y, y, &Ax_dot);
            }
            
            if (my_pid == 0) {
                printf("(  b,   b) = %.20e\n",  (double)b_dot);
                printf("(  x,   x) = %.20e\n",  (double)x_dot);
                printf("(A*b, A*b) = %.20e\n", (double)Ab_dot);
                printf("(A*x, A*x) = %.20e\n", (double)Ax_dot);
            }


            // 计算初始残差以便确定收敛阈值
            data_t r_nrm2 = 0.0, b_nrm2 = 0.0;
            // check_residual(A, x, b, y, r_nrm2, b_nrm2);
            // config.atol = r_nrm2 * 1.0e-6;

            buildup_solver(solver_name, prec_name, config);

            idx_t num_iterations = 0;
            data_t t_setup = 0.0, t_solve = 0.0, final_res_norm = 0.0;

            if (obj_type == HYPRE_SSTRUCT)
                setup_and_solve(solver_name, prec_name, A, b, x, t_setup, t_solve, final_res_norm, num_iterations);
            else if (obj_type == HYPRE_PARCSR) {
                setup_and_solve(solver_name, prec_name, par_A, par_b, par_x, t_setup, t_solve, final_res_norm, num_iterations);
                HYPRE_SStructVectorGather(b);
                HYPRE_SStructVectorGather(x);
            }

            // 计算真实残差
            check_residual(prec_name, A, x, b, y, r_nrm2, b_nrm2);
            if (my_pid == 0) {
                printf("\033[1;35mtrue ||r|| = %20.16e ||r||/||b||= %20.16e\033[0m\n", r_nrm2, r_nrm2/b_nrm2);
                printf("Iterations = %d\n", num_iterations);
                printf("Time cost %.5f %.5f %.5f %d\n", (double)t_setup, (double)t_solve, (double)(t_setup + t_solve), num_iterations);
                printf("Final Relative Residual Norm = %e\n", (double)final_res_norm);
                printf("\n");
            }

            records[test].iter = num_iterations;
            records[test].setup = t_setup;
            records[test].solve = t_solve;

            destroy_solver(solver_name, prec_name);
        }// test loop

        if (my_pid == 0) {
            std::sort(records, records + TEST_CNT);
            const TEST_RECORD & best = records[0], & worst = records[TEST_CNT - 1];
            TEST_RECORD avg;
            for (idx_t i = 0; i < TEST_CNT; i++) {
                avg.setup += records[i].setup;
                avg.solve += records[i].solve;
                avg.prec  += records[i].prec;
                avg.iter  += records[i].iter;
            }
            avg.setup /= TEST_CNT;
            avg.solve /= TEST_CNT;
            avg.prec  /= TEST_CNT;
            avg.iter  /= TEST_CNT;
            data_t min_tot = best .setup + best .solve;
            data_t max_tot = worst.setup + worst.solve;
            data_t avg_tot = avg  .setup + avg  .solve;
            printf("\n  Summary of %d tests:\n", TEST_CNT);
            printf("     Setup time  Solve time  #Iter  Total time\n");
            printf("Min    %.6f    %.6f  %5d  %.6f\n", best .setup, best .solve, best .iter, min_tot);
            printf("Avg    %.6f    %.6f  %5d  %.6f\n", avg  .setup, avg  .solve, avg  .iter, avg_tot);
            printf("Max    %.6f    %.6f  %5d  %.6f\n", worst.setup, worst.solve, worst.iter, max_tot);    
        }

        HYPRE_SStructMatrixDestroy(A);
        HYPRE_SStructVectorDestroy(b); HYPRE_SStructVectorDestroy(x); HYPRE_SStructVectorDestroy(y); HYPRE_SStructVectorDestroy(x0); 
        HYPRE_SStructGraphDestroy(ssgraph);
        HYPRE_SStructStencilDestroy(stencil);
        HYPRE_SStructGridDestroy(ssgrid);
    }
    HYPRE_Finalize();
    MPI_Finalize();
    return 0;
}